
#include <cstdlib>
#include <iostream>
#include <sstream>
#include <fstream>
#include <string>

#include "../core/tconstraints.h"
#include "../meshdef/mesh.h"
#include "../meshdef/mlsMeshDeformer.h"
#include "../meshdef/mlsMeshDeformerFS.h"
#include "../meshdef/mlsMeshDeformerVS.h"
#include "../glslKernel.h"

#include <GL/glut.h>
#include "../arcball.h"
#include "../atmatrix.h"
#include "../selection.h"
#include "../glDrawFunctions.h"
#include "../typedefs.h"

//~ #include "../mlsTimes/mlsRaw.h"
#include "../mlsTimes/mlsComparison.h"

using namespace std;

int winw = 700, winh = 700; // window size
int xini, yini; // initial mouse position
int buttonpressed; // 
MeshRenderMode rendermode; // mesh render mode
InteractionMode interactionmode; // interaction mode

Point3 grabStart, grabLast;
int id_moving;

AMatrix<GLfloat>     sceneT;
AMatrix<GLfloat>     sceneIniT;
ArcBall              arcball;
Selection            selection;

typedef TConstraints<Point3>                 PConstraints;
typedef MlsMeshDeformer<PConstraints>        MLSMeshDeformer;
Mesh                                         *mesh=0;
MLSMeshDeformer                              *deformer=0;
PConstraints                                 constraints;
PConstraints                                 constraintsIni;

MlsMeshDeformerVS *mlsdeformerVS=0;
MlsMeshDeformerFS *mlsdeformerFS=0;

vector<Point3> positionsIni;

GLuint vbo; /// vertex buffer object id
GLuint nbo; /// normal buffer object id
GLuint ibo; /// element buffer object id

/// time info
float      cfps;
time_t     start, end;
unsigned   frame;
/// Info variables
void* font=GLUT_BITMAP_TIMES_ROMAN_24;
void* fontsmall=GLUT_BITMAP_HELVETICA_12;
string nvertices("Vertices: ");
string nfaces("Faces: ");
string fps("FPS: ");
string options("Options: p/w/h/l/f/s/+/-");
string mode("Mode: Viewing");
string nconstraints("constraints: ");

/// Renders info about the model
void renderInfo();
/// Returns the translation vector in world coordinates of two points in screen space
Vector3 screenToWorld(int xi, int yi, int xe, int ye);
/// Returns the world coordinates of a point in screen space
Point3 screenToWorld(int x, int y);
/// Implements the "zoomin" operation. Narrows the view of the scene.
void zoomIn();
/// Implements the "zoomout" operation. Enlarges the view of the scene.
void zoomOut();

/// Defines the processor used for computing the deformation
char PROCESSOR = 'C'; // C: CPU, V: vertex shader, F: fragment shader
bool flagCtrlPoints = true;

/// OpenGL display function
void display(void) {
   
   glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
   glMatrixMode(GL_PROJECTION);
   glLoadIdentity();
   glOrtho(-1.0,1.0,-1.0,1.0,-4,4);
   
   //~ glOrtho(-5,5,-5,5,-20,20);
   //~ gluLookAt(1,0,0,0,0,0,0,0,1);
   
   //~ float nearplane=20, farplane=400;
   //~ gluPerspective(40*1.0, (GLfloat)winw/(GLfloat)winh, nearplane, farplane);
   //~ gluLookAt(0,-150,90,0,0,0,0,0,1);
   
   glMatrixMode(GL_MODELVIEW);
   glLoadIdentity();
   glViewport(0, 0, (GLsizei) winw, (GLsizei) winh);
   
   glMultMatrixf(&sceneT);
   
   if (interactionmode == DEFORMATION_MODE) {
      if (PROCESSOR == 'C' or PROCESSOR == 'F') {
         glDrawModel(mesh, rendermode);
      }
      else if (PROCESSOR == 'V') {
         mlsdeformerVS->use();
         glDrawModel(mesh, rendermode);
         mlsdeformerVS->use(false);
      }
   } else {
      glDrawModel(mesh, rendermode);
   }
   
   /// draw constraints
   if (flagCtrlPoints) {
      glPushAttrib(GL_COLOR_BUFFER_BIT | GL_LIGHTING_BIT | GL_POINT_BIT);
      glDisable(GL_LIGHTING);
      glEnable(GL_POINT_SMOOTH);
      glEnable(GL_BLEND);
      glColor4f(0,0,0,0);
      glPointSize(12);
      glDrawConstraints(constraints);
      glPointSize(1.0);
      glPopAttrib();
   }
   
   /// FPS info
   frame++;
   std::time(&end);
   double ddifftime = std::difftime(end, start);
   if (ddifftime > 1) {
      cfps = frame/ddifftime;
      start = end;
      frame = 0;
      ostringstream s1;
      s1<<cfps;
      fps = "FPS: " + s1.str();
   }
   
   glPushMatrix();
   glLoadIdentity();
   renderInfo();
   glPopMatrix();
   
   glutSwapBuffers();
}

/// Mouse press function, 
void mouseclick(int button, int state, int x, int y) {
   
   id_moving = -1;
   if (state == GLUT_DOWN) {
      xini=x; yini=winh-y;
      buttonpressed = button;
      if (buttonpressed == GLUT_LEFT_BUTTON and 
         interactionmode == DEFORMATION_MODE) {
         //~ cout<<"buttonpressed == GLUT_LEFT_BUTTON"<<endl;
         
         std::list<unsigned int> selected;
         selection.begin(x, winh-y);
         glDrawConstraints(constraints);
         selection.end(selected);
         if (not selected.empty()) id_moving = selected.front();
         
         grabLast = grabStart = screenToWorld(x, winh-y);
      }
      else if (buttonpressed == GLUT_MIDDLE_BUTTON) {
         arcball.click(xini,yini);
      }
      else if (buttonpressed == GLUT_RIGHT_BUTTON) {
         grabLast = grabStart = screenToWorld(x, winh-y);
      }
      sceneIniT = sceneT;
   }
   else if (state == GLUT_UP) {
      //~ cout<<"buttonpressed:GLUT_UP == GLUT_LEFT_BUTTON"<<endl;
      if (buttonpressed == GLUT_LEFT_BUTTON 
         and interactionmode == VIEW_MODE) {
         
         y = winh-y;
         std::list<unsigned int> selected;
         if (abs(xini-x)<4 or abs(yini-y)<4) {
            selection.begin(x, y);
            glDrawModelFaces(mesh);
            selection.end(selected);
            // Include the object being altered
            if (not selected.empty()) {
               unsigned idf = selected.front();
               Facet_iterator fi = mesh->facets_begin();
               advance(fi, idf);
               
               unsigned idv = fi->halfedge()->vertex()->id();
               cout<<"idv point: "<<idv<<endl;
               const Point3 p = fi->halfedge()->vertex()->point();
               constraints.push_back(p, idv);
            }
         }
         else {
            selection.begin((int)(xini+0.5*(x-xini)), (int)(yini+0.5*(y-yini)), abs(x-xini), abs(y-yini));
            glDrawModelVertices(mesh);
            selection.endAll(selected);
            
            std::list<unsigned int>::iterator it = selected.begin();
            while (it != selected.end()) {
               unsigned idv = *it++;
               Vertex_iterator vi = mesh->vertices_begin();
               advance(vi, idv);
               cout<<idv<<endl;
               constraints.push_back(vi->point(), vi->id());
            }
         }
         // label
         ostringstream s2;
         s2<<constraints.size();
         nconstraints = "constraints: " + s2.str();
      }
   }
   glutPostRedisplay();
}

AMatrix<float> mT; /// an auxiliar matrix
/// Mouse move function
void mousemove(int x, int y) {
   
   if (buttonpressed == GLUT_LEFT_BUTTON and 
      interactionmode == DEFORMATION_MODE and id_moving >= 0) {
      
      Point3 grabEnd = screenToWorld(x, winh-y);
      Vector3 t = grabEnd-grabLast;
      const Point3 &p = constraints.getConstraint(id_moving);
      constraints.setConstraint(id_moving, p + t);
      
      if (PROCESSOR == 'C') { /// CPU
         //deformer->execute(constraints);
         
         
         float p[3], pd[3]; float n[3], nd[3];
         Vertex_iterator vi = mesh->vertices_begin();
         for (int i=0; i<mesh->nvertices(); ++i) {
            const Point3 &po = positionsIni[i];
            p[0] = po[0]; p[1] = po[1]; p[2] = po[2];
            mlsEval(constraintsIni.positions, constraints.positions, p, pd, n, nd);
            //~ cout<<"xxx"<<endl;
            //~ mlsEval_OM_Horn(constraintsIni.positions, constraints.positions, p, pd);
            //~ mlsEval_Quat(constraintsIni.positions, constraints.positions, p, pd);
            //~ mlsEval_SVD(constraintsIni.positions, constraints.positions, p, pd);
            //~ mlsEval_nearly(constraintsIni.positions, constraints.positions, p, pd);
            //~ mlsEval_dualquat(constraintsIni.positions, constraints.positions, p, pd);
            vi[i].point() = Point3(pd[0],pd[1],pd[2]);
         }
         
      }
      else if (PROCESSOR == 'V') { /// GPU: vertex shader
         mlsdeformerVS->execute(constraints);
      }
      else if (PROCESSOR == 'F') { /// GPU: fragment shader
         mlsdeformerFS->execute(constraints);
      }
      
      grabLast = grabEnd;
   }
   else if (buttonpressed == GLUT_MIDDLE_BUTTON) { // rotation handler
      arcball.drag(x, winh-y, &mT);
      sceneT = mT*sceneIniT;
   }
   else if (buttonpressed == GLUT_RIGHT_BUTTON) { // translation handler
      Vector3 t = screenToWorld(xini, yini, x, winh-y);
      
      mT.identity();
      mT.translation(t[0], t[1], t[2]);
      //~ sceneT = mT*sceneIniT;
      sceneT = sceneIniT*mT;
   }
   glutPostRedisplay();
}

vector<int> handlesids; // ids of control points handles
Number MT[4][4];
float angle=0.0;

/// OpenGL key press function
void keypress(unsigned char key, int x, int y) {
   
   switch (key) {
      case 'p':
      case 'P': rendermode = POINTS; break;
      case 'w':
      case 'W': rendermode = WIREFRAME; break;
      case 'h':
      case 'H': rendermode = HIDDENLINES; break;
      case 'l':
      case 'L': rendermode = FLATLINES; break;
      case 'f':
      case 'F': rendermode = FLAT; break;
      case 's':
      case 'S': rendermode = SMOOTH; break;
      case 'c':
      case 'C': flagCtrlPoints = not flagCtrlPoints; break;
      case 'g':
      case 'G': {
         mesh->saveShapeOFF("outputMesh.off");
         break;
      }
      
      case 'd': {
         interactionmode = DEFORMATION_MODE;
         mode = "Mode: Deformation";
         angle = 0;
         
         constraintsIni = constraints;
         if (PROCESSOR == 'C') { /// CPU
            if (deformer!=0) delete deformer;
            deformer = new MLSMeshDeformer(mesh, constraints, false);
            
            Vertex_iterator vi = mesh->vertices_begin();
            for (int i=0; i<mesh->nvertices(); ++i) {
               positionsIni.push_back(vi[i].point());
            }
            
            /// cactus
            //~ const Point3 &p = constraints.getConstraint(constraints.size()-1);
            //~ Vector3 t(-1.3,-0.7,0.0);
            //~ constraints.setConstraint(constraints.size()-1, p + t);
            //~ deformer->execute(constraints);
            
            /// homer
            //~ const Point3 &p = constraints.getConstraint(constraints.size()-1);
            //~ Vector3 t(0.0,-0.8,1.2);
            //~ constraints.setConstraint(constraints.size()-1, p + t);
            //~ deformer->execute(constraints);
            
            /// dino
            //~ const Point3 &p = constraints.getConstraint(constraints.size()-1);
            //~ Vector3 t(0.4,1.5,0.0);
            //~ constraints.setConstraint(constraints.size()-1, p + t);
            //~ deformer->execute(constraints);
            
            /// girl
            //~ const Point3 &p = constraints.getConstraint(0);
            //~ Vector3 t(-0.1,-0.1,0.0);
            //~ constraints.setConstraint(0, p + t);
            //~ deformer->execute(constraints);
            
            /// comparisons
            //~ const Point3 &p = constraints.getConstraint(0);
            // Vector3 t(1.0,-1.0,0.0); // girln
            // Vector3 t(0.0,-0.8,0.0); // horse
            //~ Vector3 t(0.0,0.4,0.0); // dragon
            //~ constraints.setConstraint(0, p + t);
            /*
            const Point3 &p0 = constraints.getConstraint(0);
            const Point3 &p1 = constraints.getConstraint(1);
            Vector3 t(0.0,0.4,0.0); // dragon
            constraints.setConstraint(0, p0 + t);
            constraints.setConstraint(1, p1 - t/2.0);
            
            deformer->execute(constraints);
            */
         }
         else if (PROCESSOR == 'V') { /// GPU: vertex shader
            if (mlsdeformerVS!=0) delete mlsdeformerVS;
            mlsdeformerVS = new MlsMeshDeformerVS(constraints);
         }
         else if (PROCESSOR == 'F') { /// GPU: fragment shader
            if (mlsdeformerFS!=0) delete mlsdeformerFS;
            mlsdeformerFS = new MlsMeshDeformerFS(mesh, constraints);
         }
         
         break;
      }
      
      case 'r': {
         angle += 10.0;
         float a = angle*3.1416/180.0; // angle
         //~ MT[0][0]=cos(a);  MT[0][1]=0.0; MT[0][2]=sin(a); MT[0][3]=0.0;
         //~ MT[1][0]=0.0;     MT[1][1]=1.0; MT[1][2]=0.0;    MT[1][3]=0.0;
         //~ MT[2][0]=-sin(a); MT[2][1]=0.0; MT[2][2]=cos(a); MT[2][3]=0.0;
         //~ MT[3][0]=0.0;     MT[3][1]=0.0; MT[3][2]=0.0;    MT[3][3]=1.0;
         
         MT[0][0]=cos(a); MT[0][1]=-sin(a); MT[0][2]=0.0; MT[0][3]=0.0;
         MT[1][0]=sin(a); MT[1][1]=cos(a);  MT[1][2]=0.0; MT[1][3]=0.0;
         MT[2][0]=0.0;    MT[2][1]=0.0;     MT[2][2]=1.0; MT[2][3]=0.0;
         MT[3][0]=0.0;    MT[3][1]=0.0;     MT[3][2]=0.0; MT[3][3]=1.0;
         
         for (int i=0; i<handlesids.size(); ++i) {
            const Point3 &p = constraintsIni.getConstraint(handlesids[i]);
            Number x = p[0]*MT[0][0] + p[1]*MT[0][1] + p[2]*MT[0][2] + MT[0][3];
            Number y = p[0]*MT[1][0] + p[1]*MT[1][1] + p[2]*MT[1][2] + MT[1][3];
            Number z = p[0]*MT[2][0] + p[1]*MT[2][1] + p[2]*MT[2][2] + MT[2][3];
            constraints.setConstraint(handlesids[i], Point3(x,y,z));
         }
         deformer->execute(constraints);
         break;
      }
      
      case '+': { zoomIn(); break; }
      case '-': { zoomOut(); break; }
      case 27: {
         if (mesh!=0) delete mesh; exit(1);
      }
   }
   glutPostRedisplay();
}

/// auxiliar variables
GLdouble mMtx[16], pMtx[16];  /// modelview/projection matrix
GLint viewport[4];            /// the viewport

/// Returns the world coordinates of a point in screen space
Point3 screenToWorld(int x, int y) {
   
   glGetIntegerv(GL_VIEWPORT, viewport);
   glGetDoublev(GL_PROJECTION_MATRIX, pMtx);
   glGetDoublev(GL_MODELVIEW_MATRIX, mMtx);
   
   GLdouble xw, yw, zw;
   gluUnProject(x, y, 1.0, mMtx, pMtx, viewport, &xw, &yw, &zw);
   return Point3(xw, yw, zw);
}

/// Returns the translation vector in world coordinates of two points in screen space
Vector3 screenToWorld(int xi, int yi, int xe, int ye) {
   
   glGetIntegerv(GL_VIEWPORT, viewport);
   glGetDoublev(GL_PROJECTION_MATRIX, pMtx);
   glGetDoublev(GL_MODELVIEW_MATRIX, mMtx);
   
   GLdouble ws[3], we[3];
   gluUnProject(xi, yi, 0, mMtx, pMtx, viewport, &ws[0], &ws[1], &ws[2]);
   gluUnProject(xe, ye, 0, mMtx, pMtx, viewport, &we[0], &we[1], &we[2]);
   
   return Vector3(we[0]-ws[0], we[1]-ws[1], we[2]-ws[2]);
}

/// OpenGL reshape funcion
void reshape (int w, int h) {
   
   glMatrixMode(GL_PROJECTION);
   glLoadIdentity();
   glOrtho(-1,1,-1,1,-4,4);
   
   glMatrixMode(GL_MODELVIEW);
   glLoadIdentity();
   winw = w; winh = h;
   glViewport (0, 0, (GLsizei) w, (GLsizei) h);
   glutPostRedisplay();
   
   arcball.setBounds(winw, winh);
}

/// OpenGL initializations
void init(void) {
   
   /// //////////////////////
   glEnable(GL_NORMALIZE);
   glEnable(GL_DEPTH_TEST);
   glEnable(GL_DITHER);
   glEnable(GL_LIGHTING);
   glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, GL_FALSE);
   
   //~ GLfloat mat_A[] = {0.3, 0.3, 0.3, 1.0};
   //~ GLfloat mat_D[] = {1.0, 0.4, 0.0, 1.0};
   //~ GLfloat mat_S[] = {1.0, 1.0, 1.0, 1.0};
   //~ GLfloat shine[] = {128.0};
   GLfloat mat_A[] = {0.5, 0.5, 0.5, 1.0};
   GLfloat mat_D[] = {0.0, 0.5, 0.0, 1.0};
   GLfloat mat_S[] = {1.0, 1.0, 1.0, 1.0};
   GLfloat shine[] = {128.0};
   glMaterialfv(GL_FRONT, GL_AMBIENT,   mat_A);
   glMaterialfv(GL_FRONT, GL_DIFFUSE,   mat_D);
   glMaterialfv(GL_FRONT, GL_SPECULAR,  mat_S);
   glMaterialfv(GL_FRONT, GL_SHININESS, shine);
   
   GLfloat pos_L0[] = { 10.0, 10.0, 0.0, 0.0};
   GLfloat pos_L1[] = {-10.0, 10.0, 0.0, 0.0};
   GLfloat pos_L2[] = { 0.0, 0.0, 10.0, 0.0};
   GLfloat col_2[]  = { 1.0, 1.0, 1.0, 1.0};
   glEnable(GL_LIGHT0);
   glEnable(GL_LIGHT1);
   glEnable(GL_LIGHT2);
   glLightfv(GL_LIGHT0, GL_POSITION, pos_L0);
   glLightfv(GL_LIGHT1, GL_POSITION, pos_L1);
   glLightfv(GL_LIGHT2, GL_POSITION, pos_L2);
   glLightfv(GL_LIGHT2, GL_DIFFUSE,  col_2);
   glLightfv(GL_LIGHT2, GL_SPECULAR, col_2);
   /// 
   
   glClearColor(1,1,1,1);
   glColor4f(1,1,0,1);
   rendermode = SMOOTH;
   interactionmode = VIEW_MODE;
   
   arcball.setBounds(winw, winh);
   sceneT.identity();
   
   std::time(&start);
   frame = 0;
}

/// Implements the "zoomin" operation. Narrows the view of the scene.
void zoomIn() {
   
   AMatrix<float> mT;
   mT.identity();
   mT.scaling(1.2);
   sceneT = mT*sceneT;
}

/// Implements the "zoomout" operation. Enlarges the view of the scene.
void zoomOut() {
   
   AMatrix<float> mT;
   mT.identity();
   mT.scaling(1.0/1.2);
   sceneT = mT*sceneT;
}

/// Renders a string
void renderBitmapString(float x, float y, float z, void *font, const char *string) {
   const char *c;
   glRasterPos3f(x, y,z);
   for (c=string; *c != '\0'; c++) {
      glutBitmapCharacter(font, *c);
   }
}

/// Renders info about the model
void renderInfo() {
   
   glDisable(GL_LIGHTING); glDisable(GL_DEPTH_TEST);
   glColor4f(0,0,0,1);
   renderBitmapString(-0.9,0.90,0,(void *)fontsmall, nvertices.c_str());
   renderBitmapString(-0.9,0.85,0,(void *)fontsmall, nfaces.c_str());
   renderBitmapString(-0.9,0.80,0,(void *)fontsmall, fps.c_str());
   renderBitmapString(-0.9,0.75,0,(void *)fontsmall, options.c_str());
   renderBitmapString(-0.9,0.65,0,(void *)fontsmall, nconstraints.c_str());
   glEnable(GL_LIGHTING); glEnable(GL_DEPTH_TEST);
}

void printInfoAndExit() {
   
   cout<<endl;
   cout<<"Usage: meshdef filename.off"<<endl;
   exit(0);
}

void load_control_points() {
   
   ifstream file("cactus.sel");
   if (not file.is_open()) {
      cout<<"I/O error!!!"<<endl;
      assert(false);
   }
   int idv=0; int t=0;
   while(not file.eof()) {
      file>>t;
      if (t == 0 or t == 2) {
         Vertex_iterator vi = mesh->vertices_begin();
         advance(vi, idv);
         constraints.push_back(vi->point(), vi->id());
         if (t==2) handlesids.push_back(constraints.size()-1);
      }
      //~ else if (t == 1) {      }
      if (idv == 4800) cout<<"tipo: "<<t<<" ::: "<<constraints.size()-1<<endl;
      idv++;
   }
   /// set a control point on the top
   //~ Vertex_iterator vi = mesh->vertices_begin();
   //~ constraints.push_back(vi[4558].point(), vi[4558].id());
   
   cout<<"number of control points: "<<constraints.size()<<endl;
   
   /// load the transformation
   ifstream fileT("cactus.def");
   if (not fileT.is_open()) {
      cout<<"I/O error!!!"<<endl;
      assert(false);
   }
   fileT>>MT[0][0]>>MT[0][1]>>MT[0][2]>>MT[0][3];
   fileT>>MT[1][0]>>MT[1][1]>>MT[1][2]>>MT[1][3];
   fileT>>MT[2][0]>>MT[2][1]>>MT[2][2]>>MT[2][3];
   fileT>>MT[3][0]>>MT[3][1]>>MT[3][2]>>MT[3][3];
}

void open_mesh(string filename) {
   
   string extensao =filename.substr(filename.size()-3,filename.size());
   if (extensao == "off") {
      mesh = Mesh::load_off(filename.c_str());
      ostringstream s1, s2;
      s1<<mesh->nvertices();
      s2<<mesh->nfaces();
      nvertices += s1.str();
      nfaces += s2.str();
   }
   
   // homer
   //~ Vertex_iterator vi = mesh->vertices_begin();
   //~ constraints.push_back(vi[600].point(), vi[3089].id());
   //~ constraints.push_back(vi[4173].point(), vi[4173].id());
   //~ constraints.push_back(vi[3209].point(), vi[3209].id());
   //~ constraints.push_back(vi[3218].point(), vi[3218].id());
   //~ constraints.push_back(vi[2816].point(), vi[2816].id());
   //~ constraints.push_back(vi[2938].point(), vi[2938].id());
   //~ constraints.push_back(vi[734].point(), vi[734].id());
   //~ constraints.push_back(vi[4068].point(), vi[4068].id());
   //~ constraints.push_back(vi[2974].point(), vi[2974].id());
   //~ constraints.push_back(vi[2973].point(), vi[2973].id());
   //~ constraints.push_back(vi[3961].point(), vi[3961].id());
   
   /*
   // cactus
   Vertex_iterator vi = mesh->vertices_begin();
   constraints.push_back(vi[3089].point(), vi[3089].id());
   constraints.push_back(vi[1852].point(), vi[1852].id());
   constraints.push_back(vi[3094].point(), vi[3094].id());
   constraints.push_back(vi[5032].point(), vi[5032].id());
   constraints.push_back(vi[3864].point(), vi[3864].id());
   constraints.push_back(vi[4949].point(), vi[4949].id());
   constraints.push_back(vi[349].point(), vi[349].id());
   constraints.push_back(vi[3054].point(), vi[3054].id());
   constraints.push_back(vi[571].point(), vi[571].id());
   constraints.push_back(vi[3674].point(), vi[3674].id());
   constraints.push_back(vi[2138].point(), vi[2138].id());
   constraints.push_back(vi[4558].point(), vi[4558].id());
   */
   
   // dino
   //~ Vertex_iterator vi = mesh->vertices_begin();
   //~ constraints.push_back(vi[9862].point(), vi[9862].id());
   //~ constraints.push_back(vi[4230].point(), vi[4230].id());
   //~ constraints.push_back(vi[2221].point(), vi[2221].id());
   //~ constraints.push_back(vi[4156].point(), vi[4156].id());
   //~ constraints.push_back(vi[13780].point(), vi[13780].id());
   
   /// ///
   //~ load_control_points();
   
   /// comparison data
   // girl
   //~ int i0=8815, i1=6077, i2=2288, i3=3068, i4=7179;
   //~ Vertex_iterator vi = mesh->vertices_begin();
   //~ constraints.push_back(vi[i0].point(), vi[i0].id());
   //~ constraints.push_back(vi[i1].point(), vi[i1].id());
   //~ constraints.push_back(vi[i2].point(), vi[i2].id());
   //~ constraints.push_back(vi[i3].point(), vi[i3].id());
   //~ constraints.push_back(vi[i4].point(), vi[i4].id());
   
   // horse
   //~ int i0=15983, i1=12, i2=32, i3=220, i4=11361, i5=10461, i6=11956, i7=10756, i8=97;
   //~ Vertex_iterator vi = mesh->vertices_begin();
   //~ constraints.push_back(vi[i0].point(), vi[i0].id());
   //~ constraints.push_back(vi[i1].point(), vi[i1].id());
   //~ constraints.push_back(vi[i2].point(), vi[i2].id());
   //~ constraints.push_back(vi[i3].point(), vi[i3].id());
   //~ constraints.push_back(vi[i4].point(), vi[i4].id());
   //~ constraints.push_back(vi[i5].point(), vi[i5].id());
   //~ constraints.push_back(vi[i6].point(), vi[i6].id());
   //~ constraints.push_back(vi[i7].point(), vi[i7].id());
   //~ constraints.push_back(vi[i8].point(), vi[i8].id());
   
   // dragon
   //~ int i0=12473, i1=13693, i2=8118, i3=6911, i4=2724, i5=5968;
   //~ Vertex_iterator vi = mesh->vertices_begin();
   //~ constraints.push_back(vi[i0].point(), vi[i0].id());
   //~ constraints.push_back(vi[i1].point(), vi[i1].id());
   //~ constraints.push_back(vi[i2].point(), vi[i2].id());
   //~ constraints.push_back(vi[i3].point(), vi[i3].id());
   //~ constraints.push_back(vi[i4].point(), vi[i4].id());
   //~ constraints.push_back(vi[i5].point(), vi[i5].id());
   
   /*
   // bar model: bar0high2.off or bar0high3.off
   ifstream file("bar0high2o3_128.cp");
   int id;
   Vertex_iterator vi = mesh->vertices_begin();
   while(not file.eof()) {
      file>>id;
      constraints.push_back(vi[id].point());
   }
   // insert the non-fixed handles
   for (int i=0; i<constraints.size()/2; ++i)
      handlesids.push_back(i);
   */
   
   ostringstream s2;
   s2<<constraints.size();
   nconstraints = "constraints: " + s2.str();
   
   assert(mesh != 0);
}

/// Process the program's parameters
void input(int argc, char **argv) {
   
   int nargs = argc; // 
   if (nargs<2 or nargs>2) printInfoAndExit();
   
   open_mesh(argv[1]);
}

int main(int argc, char **argv) {
   
   // Process the input parameters
   input(argc, argv);
   
   // Glut setup
   glutInit(&argc, argv);
   glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);
   glutInitWindowSize(winw, winh);
   glutInitWindowPosition(150, 0);
   glutCreateWindow("MLS mesh deformation prototype");
   
   // OpenGL and data initialization
   init();
   
   // Callbacks registration
   glutDisplayFunc(display);
   glutReshapeFunc(reshape);
   glutMouseFunc(mouseclick);
   glutMotionFunc(mousemove);
   glutKeyboardFunc(keypress);
   glutMainLoop();
   
   return 0;
}
